Individual household electric power consumption data: Exploratory Analysis¶
Ripunjoy Gohain
ripunjoygohain79@gmail.com
Data Information:
- date: Date in format dd/mm/yyyy
- time: time in format hh:mm:ss
- global_active_power: household global minute-averaged active power (in kilowatt)
- global_reactive_power: household global minute-averaged reactive power (in kilowatt)
- voltage: minute-averaged voltage (in volt)
- global_intensity: household global minute-averaged current intensity (in ampere)
- sub_metering_1: energy sub-metering No. 1 (in watt-hour of active energy). It corresponds to the kitchen, containing mainly a dishwasher, an oven and a microwave (hot plates are not electric but gas powered).
- sub_metering_2: energy sub-metering No. 2 (in watt-hour of active energy). It corresponds to the laundry room, containing a washing-machine, a tumble-drier, a refrigerator and a light.
- sub_metering_3: energy sub-metering No. 3 (in watt-hour of active energy). It corresponds to an electric water-heater and an air-conditioner.
In [1]:
# coding: utf-8
# !/usr/bin/env python3
import sys
import pandas as pd
pd.set_option("display.max_columns", None)
pd.set_option('float_format', '{:4f}'.format)
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
from scipy import stats
In [2]:
print("Python Version: ", sys.version)
print("Pandas Version: ", pd.__version__)
print("Numpy Version: ", np.__version__)
print("Seaborn Version: ", sns.__version__)
In [3]:
# dataset is txt form, deliminated with ";", reading the dataset variables as string so setting the low_memory to false
df = pd.read_csv("./data/household_power_consumption.zip", compression="zip", sep=";", low_memory=False, header=0,
infer_datetime_format=True, parse_dates={"local_time": [0, 1]}, index_col=["local_time"])
In [4]:
# converting the column headers to small letters
df.columns = [x.lower() for x in df.columns]
# missing values are represented as "?"
df.replace(to_replace="?", value=np.nan, inplace=True)
# convert all values to float32
df = df.astype(np.float32)
In [5]:
df.head(10)
Out[5]:
In [6]:
print("Shape of the data frame: ", df.shape)
In [7]:
# Calculating expected number of index for 1 minutes frequency in the provided time period
print("Expected number of index(data) in the time range: {:.0f}".format((df.index.max()-df.index.min()).total_seconds() / (60) + 1))
Shape and expected index matching, so there can be missing values but no missing index¶
In [8]:
print("Data start time: ", df.index.min())
print("Data end time: ", df.index.max())
print("Total time period: ", df.index.max()-df.index.min())
In [9]:
df.info()
In [10]:
df.describe()
Out[10]:
In [11]:
# Missing value statistics
print("Column wise missing values: ")
df.isna().sum()
Out[11]:
In [12]:
# If we drop row if all of the values are missing then remaining shape will be: Just calculating, not assigning to the dataframe
df.dropna(axis=0, how="all").shape
Out[12]:
In [13]:
# Let's check if this difference of dropping any is matching indidividual missing value count or not
df.shape[0] - df.dropna(axis=0, how="any").shape[0]
Out[13]:
In [14]:
# From above, as it's matching we can say that, if there if one tag is missing, then all other tags are missing, missing value
# is may not from some specific sensor but feels like communication erro5
print("Missing values percentage: {:.2f} %".format((df.shape[0] - df.dropna(axis=0, how="any").shape[0]) * 100 / df.shape[0]))
In [15]:
# let's check if missing values are from random timestamps or it occurs for longer time interval
missing_stat = df.dropna().index.to_series().diff().dropna().dt.total_seconds().div(60)
In [16]:
# Here actually we are calculating time difference, so diff 1 minutes means no missing value
# 2 means, 1 missing value before some particular time stamp, slly 10 means 9 missing values before that particular time
missing_stat.describe()
Out[16]:
In [17]:
print("Count of missing timestamps instances: ", missing_stat[missing_stat >= 2].shape)
missing_stat[missing_stat >= 2].sort_values(ascending=False).head(10)
Out[17]:
In [18]:
missing_stat[missing_stat >= 2].plot(figsize=(10,5), title="Missing values time range");
In [19]:
print("Consecutive missing values can go upto: {:.2f} days".format((missing_stat.max() - 1) / (60*24)))
Ways to fill missing values¶
- When missing values are coming for few minutes (say 10, 20 minutes), Linearly interpolating the values will make sense as we don't expect very drastic change in over few minutes.
- When it is in hours, we might need to look for previous few days that particular time as well as nearby times can come up with an strategy.
- In our case it is going for days (5 days). It is little difficult to fill those values, anything can happen in those 5 days. Strategies can be followed such:
- Looks for that particular weekday and time in last 30-60 days and fill the values with aggregated values from lookup time horizon.
- As we have multiple years data, look for those particular date time over years (might not work much well, year to year this can change, when houshold buys new devices expectations are that consumptions will change).
- We can't apply methods like filling missing values by looking on the other values and trying to fit some models (say knn, regression) here. Because when 1 parameter is missing, other parameters are also missing.
Due to time constraints, I am not going to experiment with these methods. We need to create a fuzzy logic based on the period of missing values To keep it simple, I will try to just fill the values by looking previous days same time power consumption. For longer time periods (like 5 days), It will just repeat for all the days. Basically the data will be filled up with last available time stamps data
In [20]:
def fill_missing_values_simple(dataFrame):
"""
Input: dataFrame
Simple approach to fill missing values
If today some value is missing, look for yesterdays same time data and fill it up with that
Output: dataFrame
"""
one_day_minutes = 60 * 24
for row in range(dataFrame.values.shape[0]):
for col in range(dataFrame.values.shape[1]):
if np.isnan(dataFrame.values[row, col]):
dataFrame.values[row, col] = dataFrame.values[row - one_day_minutes, col]
return dataFrame
In [21]:
df = fill_missing_values_simple(dataFrame=df)
In [22]:
df.shape
Out[22]:
In [23]:
# Verify if missing values are being filled
df.isna().sum()
Out[23]:
In [24]:
df.head()
Out[24]:
In [25]:
[x for x in df.columns if "metering" in x]
Out[25]:
In [26]:
# global_active_power is in kw, sub metering is in wH
# We will get other part energy being consumed by converting global_active_power to energy and subtracting all other sub metering energy
df['sub_metering_other'] = (df.global_active_power * 1000 / 60) - (df[[x for x in df.columns if "metering" in x]].sum(axis=1))
In [27]:
df.sub_metering_other.describe()
Out[27]:
In [28]:
df[df.sub_metering_other<0].shape
Out[28]:
Above observation is an interesting observation¶
We expected that global_active_power (converted to energy) should be equal or lesser than all 3 energy meters combined readings.
But here the negative values indicating that which is not true. And this is very unlikely. Though 1K records out of 2 million records is very less, we may be able to treat this as outliers or further extend the analysis to understand this. It may be because of fauly sensor readings.
For the time being, we will change this negative value to ZERO
If time permits, we will try to understand this particular scenario in detail.
In [29]:
df.loc[df.sub_metering_other<0, "sub_metering_other"] = 0
In [30]:
df.describe()
Out[30]:
In [31]:
# Saving the dataframe for later use
df.to_csv("./data/cleaned_household_power_consumption.csv")
Plotting¶
In [32]:
# Getting
df["year"] = df.index.year
df["quarter"] = df.index.quarter
df["month"] = df.index.month
df["day"] = df.index.day
# Monday 0, Sunday 6
df["weekday"]=df.index.weekday
df["weekday"] = (df["weekday"] < 5).astype(int)
Univariate distribution of observations¶
In [33]:
cols = ["global_active_power", "global_reactive_power", "voltage", "global_intensity", "sub_metering_1", "sub_metering_2",
"sub_metering_3", "sub_metering_other"]
In [34]:
for col in cols:
sns.distplot(df[col])
plt.tight_layout()
plt.savefig("./plots/" + "distribution_" + col + ".png")
plt.show()
Normality check of the data¶
KURTOSIS: describes heaviness of the tails of data distribution
Normal Distribution has a kurtosis of close to 0. If the kurtosis is greater than zero, then distribution has heavier tails. If the kurtosis is less than zero, then the distribution is light tails.
SKEWNESS: describes the symmetry of the data distribution
If the skewness is between -0.5 and 0.5, the data are fairly symmetrical. If the skewness is between -1 and – 0.5 or between 0.5 and 1, the data are moderately skewed. If the skewness is less than -1 or greater than 1, the data are highly skewed.
In [35]:
for col in cols:
print("\n*** Normality Test for: {} ***".format(col))
print("--"*30)
stat, p = stats.normaltest(df[col])
print("Statistics(z-score for kurtosis & skewness): {:.3f}, pvalue: {:.3f}".format(stat, p))
alpha = 0.05
if p > alpha:
print("{} data looks Normal: (fail to reject H0)".format(col))
else:
print("{} data doesn't looks Normal: (reject H0)".format(col))
print( 'Kurtosis of normal distribution: {}'.format(stats.kurtosis(df[col])))
print( 'Skewness of normal distribution: {}'.format(stats.skew(df[col])))
Time Series Plots:¶
In the next Notebook interactive plots being created using bokeh so that we can drill down(Note: In GitHub we will get the source code, but not the plots. It needs a server or needs to run in Local)
Minutely Plots¶
In [36]:
for col in cols:
df[col].plot(title=col + " (minute)", figsize=(15,5))
plt.tight_layout()
plt.savefig("./plots/" + "minutely_" + col + ".png")
plt.show()
Hourly Plots¶
In [37]:
# hourly plots
for col in cols:
df[col].resample("H").mean().plot(title=col + " (hour)", figsize=(15,5))
plt.tight_layout()
plt.savefig("./plots/" + "hourly_" + col + ".png")
plt.show()
Daily Plots for different parameters¶
In [38]:
# Daily plots
for col in cols:
df[col].resample("D").mean().plot(title=col + " (day)", figsize=(15,5))
plt.tight_layout()
plt.savefig("./plots/" + "daily_" + col + ".png")
plt.show()
Weekly Plots¶
In [39]:
for col in cols:
df[col].resample("W").mean().plot(title=col + " (week)", figsize=(15,5))
plt.tight_layout()
plt.savefig("./plots/" + "weekly_" + col + ".png")
plt.show()
Monthly Plots for different parameters¶
In [40]:
# For better view lets look at Monthly plots
for col in cols:
df[col].resample("M").mean().plot(title=col + " (month)", figsize=(15,5))
plt.tight_layout()
plt.savefig("./plots/" + "monthly_" + col + ".png")
plt.show()
Yearwise boxplots with respect to weekday/weekends, months, quarter¶
In [41]:
dic={0: "Weekend", 1: "Weekday"}
df["weekday_string"] = df.weekday.map(dic)
for col in cols:
plt.figure(figsize=(10,5))
sns.boxplot("year", col, hue="weekday_string", data=df, width=0.8, fliersize=3)
plt.tight_layout()
plt.savefig("./plots/" + "box_weekday_" + col + ".png")
plt.show()
2016 has less data, show medians are going high (Only december months data)¶
In [42]:
for col in cols:
plt.figure(figsize=(15,5))
sns.boxplot("year", col, hue="month", data=df, width=0.8, fliersize=3)
plt.tight_layout()
plt.savefig("./plots/" + "box_month_" + col + ".png")
plt.show()
In [43]:
for col in cols:
plt.figure(figsize=(10,5))
sns.boxplot("year", col, hue="quarter", data=df, width=0.8, fliersize=3)
plt.tight_layout()
plt.savefig("./plots/" + "box_quarter_" + col + ".png")
plt.show()
Quarter 1 & 3 showing interesting results, as well as it is notable to see the consumption box plots. Submeter 3 is for water heater and air conditioner. We can make intuitions from that.¶
Factor plots to see the trend over year¶
In [44]:
for col in cols:
plt.figure(figsize=(10,5))
sns.factorplot("year", col, hue="weekday_string", data=df, fliersize=3)
plt.tight_layout()
plt.savefig("./plots/" + "factor_weekday_" + col + ".png")
plt.show()
In [45]:
for col in cols:
plt.figure(figsize=(10,5))
sns.factorplot("year", col, hue="quarter", data=df, fliersize=3)
plt.tight_layout()
plt.savefig("./plots/" + "factor_quarter_" + col + ".png")
plt.show()
In [46]:
for col in cols:
plt.figure(figsize=(10,5))
sns.factorplot("month", col, hue="year", data=df, fliersize=3)
plt.tight_layout()
plt.savefig("./plots/" + "factor_month_year_" + col + ".png")
plt.show()
Let's check correlation among data¶
In [47]:
# minute wise data points
plt.figure(figsize=(10,8))
sns.heatmap(df[cols].corr(method="pearson"), vmin=-1, vmax=1, center=0, annot=True, fmt=".2f", linewidths=0.2);
plt.savefig("./plots/corr_heatmap_minute.png")
In [48]:
# hour wise data points
plt.figure(figsize=(10,8))
sns.heatmap(df[cols].resample("H").mean().corr(method="pearson"), vmin=-1, vmax=1, center=0, annot=True, fmt=".2f", linewidths=0.2);
plt.savefig("./plots/corr_heatmap_hour.png")
In [49]:
# day wise data points
plt.figure(figsize=(10,8))
sns.heatmap(df[cols].resample("D").mean().corr(method="pearson"), vmin=-1, vmax=1, center=0, annot=True, fmt=".2f", linewidths=0.2);
plt.savefig("./plots/corr_heatmap_day.png")
In [50]:
# month wise data points
plt.figure(figsize=(10,8))
sns.heatmap(df[cols].resample("M").mean().corr(method="pearson"), vmin=-1, vmax=1, center=0, annot=True, fmt=".2f", linewidths=0.2);
plt.savefig("./plots/corr_heatmap_day.png")
